ADA Programming Guide/ADA Programming Guide.iso/ada_gwu/adaed/math/primfunc.adb

Text File | 1996-01-30 | 9.4 KB | 335 lines

--  -----------------------------------------------------------------------
--  Title:       primitive_functions
--  Last Mod:    Thu Apr 19 11:25:43 1990
--  Author:      Vincent Broman <[e-mail]>
--   Copyright 1990 Vincent Broman
--     Permission granted to copy, modify, or compile this software for
--   one's own use, provided that this copyright notice is preserved intact.
--     Permission granted to distribute compiled binary copies of this
--   software which are linked in with some other application.
--     Permission granted to distribute other copies of this software,
--   provided that (1) any copy which is not source code, i.e. not in the
--   form in which the software is usually maintained, must be accompanied
--   by a copy of the source code from which it was compiled, and (2) the
--   one distributing it must refrain from imposing on the recipient
--   further restrictions on the distribution of this software.
--
--  Visibility:
--  Description:
--               functions on floating point types whose implementation
--               depends on access to the mantissa/exponent representation
--               of the floating point number.  this includes
--               integer/fraction operations.
--
--               This is a portable, slow implementation.
--
--  Exceptions:  numeric_error upon overflow of function scale,
--               constraint_error only if type Float is constrained.
--  -----------------------------------------------------------------------
package body primitive_functions is
--


   subtype expbits is integer range 0 .. 6; -- 2**(2**7) might overflow

   Float_exp: array( expbits) of integer :=
      (1, 2, 4, 8, 16, 32, 64);
   Float_power: array( expbits) of Float :=
      (2.0 **  1, 2.0 **  2, 2.0 **  4, 2.0 ** 8,
       2.0 ** 16, 2.0 ** 32, 2.0 ** 64);
   Float_neg_power: array( expbits) of Float :=
      (0.5 **  1, 0.5 **  2, 0.5 **  4, 0.5 ** 8,
       0.5 ** 16, 0.5 ** 32, 0.5 ** 64);

  -- double stuff deleted here

   function exponent( x: Float) return integer is
--
-- return the exponent k such that 1/2 <= x/(2**k) < 1,
-- or zero for x = 0.0 .
--
   begin
      if x = 0.0 then
     return 0;
      end if;
      -- nonzero x
      -- essentially, just multiply by powers of two to get in range

     declare
        ax: Float := abs(  x);
        exp: integer := 0;
        -- ax*2.0**exp is invariant
     begin
        while ax >= Float_power( expbits'last) loop
           ax := ax * Float_neg_power( expbits'last);
           exp := exp + Float_exp( expbits'last);
        end loop;
        -- ax < 2**64
        for n in reverse expbits'first .. expbits'last - 1 loop
           if ax >= Float_power( n) then
          ax := ax * Float_neg_power( n);
          exp := exp + Float_exp( n);
           end if;
           -- ax < Float_power( n)
        end loop;
        -- ax < 2
        while ax < Float_neg_power( expbits'last) loop
           ax := ax * Float_power( expbits'last);
           exp := exp - Float_exp( expbits'last);
        end loop;
        -- 2**-64 <= ax < 2
        for n in reverse expbits'first .. expbits'last - 1 loop
           if ax < Float_neg_power( n) then
          ax := ax * Float_power( n);
          exp := exp - Float_exp( n);
           end if;
           -- Float_neg_power( n) <= ax < 2
        end loop;
        -- 1/2 <= ax < 2
        if ax < 1.0 then
           return exp;
        else
           return exp + 1;
        end if;
     end;

   end exponent;


   function mantissa (x: Float) return Float is
--
-- return scale( x, - exponent( x)) if x is nonzero,  0.0 otherwise.
--
   begin
      if x = 0.0 then
     return 0.0;
      end if;
      -- nonzero x
      -- essentially, just multiply by powers of two to get in range

     declare
        ax: Float := abs( x);
     begin
        while ax >= Float_power( expbits'last) loop
           ax := ax * Float_neg_power( expbits'last);
        end loop;
        for n in reverse expbits'first .. expbits'last - 1 loop
           if ax >= Float_power( n) then
          ax := ax * Float_neg_power( n);
           end if;
        end loop;
        while ax < Float_neg_power( expbits'last) loop
           ax := ax * Float_power( expbits'last);
        end loop;
        for n in reverse expbits'first .. expbits'last - 1 loop
           if ax < Float_neg_power( n) then
          ax := ax * Float_power( n);
           end if;
        end loop;
        if ax >= 1.0 then
           ax := ax * 0.5;          -- 1 <= ax < 2
        end if;
        if x > 0.0 then
           return Float( ax);
        else
           return Float( - ax);
        end if;
     end;

   end mantissa;


   function scale (x: Float;
           k: integer) return Float is
--
-- return x * 2**k quickly, or quietly underflow to zero,
-- or raise an exception on overflow.
--
   begin
      if x = 0.0 or k = 0 then
     return x;
      end if;
      -- nonzero x and k.
      -- essentially, just multiply repeatedly by 2**(+-2**n).
      -- if is_
     declare
        y: Float  :=  x;
        exp: integer := k;
        -- y*2.0**exp is invariant
     begin
        while exp <= - Float_exp( expbits'last) loop
           y := y * Float_neg_power( expbits'last);
           exp := exp + Float_exp( expbits'last);
        end loop;
        for n in reverse expbits'first .. expbits'last - 1 loop
           if exp <= - Float_exp( n) then
          y := y * Float_neg_power( n);
          exp := exp + Float_exp( n);
           end if;
        end loop;
        -- exp >= 0
        while exp >= Float_exp( expbits'last) loop
           y := y * Float_power( expbits'last);
           exp := exp - Float_exp( expbits'last);
        end loop;
        for n in reverse expbits'first .. expbits'last - 1 loop
           if exp >= Float_exp( n) then
          y := y * Float_power( n);
          exp := exp - Float_exp( n);
           end if;
        end loop;
        -- exp = 0
        return Float( y);
     end;

   end scale;

--
-- truncation requires making all fractional bits zero in a signed
-- magnitude representation of floating point numbers.
-- for large numbers, this has no effect, as all large floats are integers.
-- for abs( x) < 1.0 this produces zero.
--

   function truncate (x: Float) return Float is
--
-- truncate x to the nearest integer value with absolute value
-- not exceeding abs( x).  No conversion to an integer type
-- is expected, so truncate cannot overflow for large arguments.
--
--
      large: Float  := 1073741824.0;
      type long is range - 1073741824 .. 1073741824;
      -- 2**30 is longer than any single-precision mantissa
      rd: Float;

   begin
      if abs( x) >= large then
     return x;
      else
     rd := Float ( long( x));
     if x >= 0.0 then
        if rd <= x then
           return rd;
        else
           return rd - 1.0;
        end if;
     else
        if rd >= x then
           return rd;
        else
           return rd + 1.0;
        end if;
     end if;
      end if;
   end truncate;


   function shorten (x: Float;
             k: integer) return Float is
--
-- set all but the k most significant bits in the mantissa of x to zero,
-- i.e. reduce the precision to k bits, truncating, not rounding.
-- shorten( x, k) = 0.0 if k < 1 and
-- shorten( x, k) = x   if k >= Float'machine_mantissa.
--
      kex: constant integer := k - exponent( x);

   begin                                -- shorten
      -- the tests on k are needed only to avoid under/overflow
      if x = 0.0 or k < 1 then
     return 0.0;
      elsif k >= Float'machine_mantissa then
     return x;
      else
     return scale( truncate( scale( x, kex)), - kex);
      end if;
   end shorten;


   function odd (x: Float) return boolean is
--
-- predicate indicates whether or not truncate( x) is an odd integer.
--
     x2: constant Float := truncate(x) / 2.0;

   begin
      return x2 /= truncate( x2);
   end odd;

--
-- all the other integer/fraction functions are based on truncate.
--

   function floor (x: Float) return Float is
--
-- return as a Float the greatest integer value <= x.
--
        sx: constant Float  := x;
        tx: constant Float := truncate( sx); -- no constraint_error
     begin
        if sx >= 0.0 or sx = tx then
           return Float( tx);
        else
           return Float( tx - 1.0); -- exactly
        end if;


   end floor;


   function ceiling (x: Float) return Float is
--
-- return as a Float the least integer value >= x.
--
        sx: constant Float :=  x;
        tx: constant Float := truncate( sx); -- no constraint_error
     begin
        if sx <= 0.0 or sx = tx then
           return Float( tx);
        else
           return Float( tx + 1.0); -- exactly
        end if;


   end ceiling;


   function round (x: Float) return Float is
--
-- return as a Float the integer value nearest x.
-- in case of a tie, prefer the even value.
--


        sx: Float  :=  x;
        tx: Float := truncate( sx);
        diff: Float := sx - tx;     -- computed exactly
     begin
        if abs( diff) < 0.5 then
           return Float( tx);
        elsif diff > 0.5 then
           return Float( tx + 1.0);
        elsif diff < -0.5 then
           return Float( tx - 1.0);
        elsif diff = 0.5 then
           if odd( tx) then
          return Float( tx + 1.0);
           else
          return Float( tx);
           end if;
        else                        -- diff = -0.5
           if odd( tx) then
          return Float( tx - 1.0);
           else
          return Float( tx);
           end if;
        end if;


   end round;

end primitive_functions;

-- $Header: p_primitive_functions_b.ada,v 3.13 90/04/19 12:34:22 broman Rel $